A Particle Method for a Collisionless Plasma with Infinite 

Mass 



Stephen Pankavich 



Department of Mathematics 
University of Texas at Arlington 
Arlington, TX 76016 
sdp@uta.edu 



Abstract 



The one-dimensional Vlasov-Poisson system is considered and a particle method is 
developed to approximate solutions without compact support which tend to a fixed 
background of charge as \x\ — > oo. Such a system of equations can be used to model 
kinetic phenomena occurring in plasma physics, such as the solar wind. The particle 
method is constructed, implemented, and used to determine information regarding 
the time asymptotics of the electrostatic field. 

Keywords: plasma, particle method, Vlasov-Poisson, infinite mass 



Introduction 

The motion of a collisionless plasma - an ionized gas of high-temperature and low- 
density - is described by the Vlasov-Maxwell system. If this fundamental kinetic 
model is posed in a two-dimensional phase space (one for spatial variables and an- 
other representing momentum) then Maxwell's equations simplify greatly and the 
problem reduces to the one-dimensional Vlasov-Poisson system. We consider this 
system of equations for the motion of negative charges upon prescribing a fixed, 
spatially-homogeneous background of positive ions: 



Here, t G [0,T] denotes time, i 6 1 is space, v e M is momentum, f(t,x,v) repre- 
sents the density of negative ions, F(v) ^ is a given function describing the fixed 
background of positive charge, and E(t,x) represents the electric field generated by 
the charges. Further define 
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to be the density of charge in the system. Unlike many classical investigations of the 
Vlasov-Poisson system, we are interested in studying how the negative ions move to 
balance the fixed positive background. Hence, we seek solutions / which tend to F 
as |x| — > oo, rather than to zero. The introduction of the fixed background implies 
that the total positive charge, total negative charge, and total energy are all infinite. 
In the present study, we will also assume neutrality 



This condition then yields zero data at infinity for the electric field, and hence a 
representation for the field results: 



To set the context of the problem, we impose assumptions on the data. Specifically, 
we will assume throughout that fo G C X (1R 2 ) is nonnegative with compact support in 
v, and there is R > such that for |x| > R, 



where F G C^(1R) is a nonnegative, even, and decreasing function. 

As one can see, no assumption of spatial compact support is made on the initial 
particle distribution f . The impetus for investigating a problem like ([T| arises from 
plasma physicists' attempts to study the stability of a two-species neutral plasma in 
which a perturbation of the distribution of negative ions is introduced from equilib- 
rium and evolves in such a way as to cancel the effects of the prescribed Maxwellian 
density of positive charge. Though our assumption of compact velocity support pre- 
cludes Maxwellian backgrounds F, the particle simulation would need to be truncated 
regardless, and the method will still approximate values of a Maxwellian. As the 
initial particle density lacks spatial decay, most of the known mathematical theory 
regarding existence and uniqueness of solutions does not apply. Recently, the author 
[7] has shown the local-in-time existence and uniqueness of solutions to ([I]) under 
similar hypothesis, and this argument can easily be adapted to extend the result to 
the assumptions above. Having answered the question of local well-posedness, a next 
logical step is to construct numerical methods that can determine the behavior of 
solutions. As such, this is the main objective of the current work. 

For problems which use a kinetic description of plasma, particle methods are often 
utilized to numerically approximate solutions and tend to be much more efficient 
and accurate than other traditional approximation techniques for partial differential 
equations, such as finite difference or finite element methods. While particle methods 
for the Vlasov-Poisson system have been previously studied [1-6], the introduction of 
non-zero spatial behavior of the particle density as \x\ — > oo eliminates the crucial 
feature of compact spatial support. Hence, any truncation of the spatial domain must 
consider the effects of particles originating from outside this region. In what follows, 





(2) 



f {x,v) = F(v) 



(3) 
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we derive and construct a method which combats this problem, while discussing both 
its abilities and limitations. 

The general structure of a particle method can be described quite simply (see pQ). 
As in other numerical methods, phase space (which is (x, v) for kinetic equations) 
is discretized into grids of finite length. Particles are then initialized with starting 
positions and velocities at time t — 0. The charge (and if necessary, current) den- 
sity is calculated from these particles, and the electric (and if necessary, magnetic) 
field is calculated from the density. Finally, we calculate the force exerted by the 
field and "move" the particles by changing their respective positions and velocities 
accordingly. Since the trajectories have now been calculated for the next time step, 
the process repeats until a stopping time t — T is reached. Weighting schemes play 
a large role in these calculations, specifically because particle charges, positions, and 
velocities must be recorded "at the particles", whereas densities, fields, and forces 
are indexed by prescribed gridpoints, with the number of particles and gridpoints 
differing dramatically. 

Typically, a particle method tracks the evolution of a finite amount of charge 
within fixed (or adaptable) spatial and velocity domains. However, for any simulation 
of 0, which models a plasma density without compact spatial support, the spatial 
domain must be truncated. Thus, we choose L > and perform the computations on 
the interval —L < x < L. The assumption on the data is reformulated to lie within 
the initial spatial domain, i.e. there is R G [0, L) such that for R < \x\ < L, 

fo(x,v) = F(v). (4) 

This ensure that the charges cancel outside of a spatial interval [-R, R] . Similarly, 
the velocity domain must be truncated. We choose Q > so that particles in the 
simulation may take on velocities only in the interval — Q < v < Q at the initial 
time. As the process continues, however, the velocity domain is enlarged to allow 
the particles to move as dictated by their interaction with the self-consistent electric 
field. The compact velocity support of / and F ensure that the largest particle 
velocity attained at any timestep is finite. Hence, at every timestep this value is 
computed and used to extend the spatial boundary of particle dynamics. Thus, the 
spatial domain is also enlarged with time, dependent upon current velocities of the 
particles, to ensure that none can escape. 

In the present context, particles are allowed to move neither into nor out of the 
computational domain. Though the positive and negative charges cancel outside 
of the spatial interval [— R, R], both positive and negative particles exist outside of 
this domain and may influence the computations. Hence, in addition to enlarging the 
spatial grid, the domain of validity, in which the particles beginning inside the interval 
[-L, L] could not have been influenced by those which began outside, is computed 
after the simulation is complete. The observed values of p and E are only considered 
valid inside this region. 
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Description of the Method 

In this section the particle method is constructed. A leap-frog scheme is utilized for 
the particle trajectories and first-order averaging methods are used to interpolate the 
field and charge density values. Begin by choosing Ax, Av > and define for every 
i,j e Z, 

Xij(0) = i- Ax 

V tJ (0)=j-Av 

q ij = f (X ij (p),V ij (p)).Ax Av. 

These quantities represent the initial particle positions, initial particle velocities, and 
the total negative charge included in the simulation, respectively. The functions Xij(t) 
and Vij(t) will be defined later for t > 0. Once they are known the approximation of 
the continuous number density is then given by 

f(t,x,v) = J2^A X ~ X v(t)) - V^t)) (5) 
where S is the Dirac mass and 5 is the first-order weighting function defined by 

5(x) = 










, if \x\ 





0, otherwise. 

Choose At > and define t n = n ■ At and xi = I ■ Ax for n E N, / G Z. We will write 

E? = E(t n , Xl ) 

for field values at the n th timestep and I th gridpoint and define the function E n (x) by 
linear interpolation of the gridpoint values E™ . 

To initiate the leap-frog scheme, we first shift the particle velocities backward by 
a half timestep using the initial field values. Hence, let 

^•r 1/2 ) = ^(o)-E(x J ,(o)).^. 

Now, for n E N U {0}, assuming Vij(t n ~ l l 2 ) has been computed, define 

s^- 1 ' 2 : = sup (v^r- 1 / 2 )! 



to be the largest velocity at time t n_1 / 2 . The length of the spatial domain at time 
zero, denoted by L, must be enlarged at every time step, so define L° = L and for 
every n E N, 

L n ._ L n-1 + gn-l/2 . Af 

Thus, L n will be the length of the spatial grid at the n th time step. This will enable 
the program to continually account for particles which began inside the inital spatial 
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domain, but due to an increase in their velocity would normally move outside of this 
region. By enlarging the spatial domain at every time step in this fashion, we are 
able to track such particles and their resulting effects on the induced field. 

Assume for n e N, Xij(t n ) and V^(£ n_1//2 ) are known for all i,j E Z, and that 

is known for \l\ < -^ L . For \l\ > in light of the neutrality of the plasma, 
we find E™ = 0. After a linear interpolation of the field values at gridpoints, define 
^■(r +1/2 ) by 




= q l] -E n {X l] (t n )). 



and 

x^(r +1 ) = x l3 (t n ) + At ■ Ki(t" +1/2 )- 

Next, we describe how to advance the electric field. Define for / e Z and iigN, 

P? +1 = J {F(v)-f(t n+ \x h v)) dv, 

and by linear interpolation of , define p n+1 (x). Then, define 

E? +1 = j X y + \y) dy. 

The process may then continue by advancing the particle positions and velocities, 
and Vij, to the next time step, t n+1 . 

Finally, as the process continues, we must determine the region on which our 
approximations are valid. Since it was necessary to truncate the spatial domain, we 
were forced to neglect the effects of particles which begin outside of our domain and 
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move with large enough velocity to enter it at some time step. Therefore, we must 
discount the approximations within the region in which these particles could have 
entered and affected the computations. Define, for every neN, 



pn ._ ^ gfc+1/2 



fc=0 



to be the total sum of largest particle velocities up to time t n+1 ^ 2 . Since no particle 
beginning outside the original spatial domain could move with velocity greater than 
P n , we may conclude that the portion of the original spatial domain that is unaffected 
by such particles at any time step t n lies inside the interval [-L + P n At, L — P n At}. 
As Figure [T] illustrates, this provides us with a specific space-time domain on which 
our numerical approximations are valid. 



Validation and Steady States 

In addition to the situation in which fo(x, v) = F(v) and hence E = 0, the previously 
described particle method was tested using a known steady state solution defined as 
follows. Let 

f(x,v)=F^-\v\ 2 + U(x)^, 

where 

_, v / -e if e < 
^ (e) = ( if e > 

and 

U(x) :=-i(l-, 2 ) 3 XH ,i)(4 
From the steady potential U(x), the resulting time- independent field is calculated by 

£ (x) := U'{x) = 3x(l — x 2 ) 2 X(-i,i)(x). 
Additionally, the charge density 

p(x) = £'(x) = 3(1 - x 2 )(l - 5a; 2 )x(-l, l)(x). 
Therefore, the distribution of positive charge is determined by p and / as 

F{x,v) = (3(1 - x 2 )(l - 5x 2 ) + ^(1 - x 2 )^j X (-iA)(x)5(v). 

Using these functions, the method was implemented for several choices of Ax, Av, 
and At in order to assure convergence to the correct steady state solution. Table 
[l] summarizes the results of these runs, listing the error found by calculating the 
difference between the known steady field solution and the computed electric field at 
every time step. 

We expect that as the mesh is refined (and the values of At, Ax, and Av decrease), 
the values of the error should also decrease at a suitable rate for each time. Choosing 
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mesh values \ time values 


t = 


t = 0.12 


t = 0.24 


t = 0.36 


t = 0.48 


At = Ax = Av = 0.04 


8.0 x lO" 3 


8.0 x 10~ 3 


0.012 


0.016 


0.018 


At = Ax = Av = 0.02 


2.0 x 10~ 3 


2.0 x 10~ 3 


3.0 x lO" 3 


0.004 


0.0005 


At = Ax = Av = 0.01 


5 x 10~ 4 


5 x 10~ 4 


8 x 10^ 4 


0.001 


0.002 


supi \E{t,xi)\ 


0.8548 


0.8568 


0.8598 


0.8620 


0.8629 



Table 1: Error of the field for steady state solution 




5 10 15 20 25 30 



Figure 2: The computed electric field (left, lined), asymptotic behavior (left, dashed), and net energy 
(right) for < t < 30 

a time t in the table and evaluating each of the three error values, we see that this 
is the case, and that as the values of spacings are halved, the error decreases by a 
factor of 4. Thus, we trust the method converges at a quadratic rate and is second 
order accurate. 

Simulation and Time Asymptotics 

Now, consider the following choices for F and /q. Let U : M 2 x [0, oo) — > R be 
defined by 

U(z,A,B) := A(B-z 2 f X{ _^ m (z), 

and then define 

F(v) = U(v,l,l) 

and 

fo(x, v) = F(v) + x ■ U(x, 1,1)- U(v, 0.1, 0.6). 
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t = 15.2 


t = 19.0 


t = 23.2 


t = 26.0 


sup \E(t,xi)\ 
i 


6.5 x 1(T 6 


5.7 x 10~ 6 


5.3 x 10~ 6 


5.0 x 10~ 6 


sup \E(t,xi) \ ■ t 
i 


9.9 x 10~ 5 


1.1 x 10" 4 


1.2 x 10" 4 


1.3 x 10" 4 



Table 2: Time rate of decay for the computed electric field 



The previously described method was implemented with this data for several 
choices of T, L, Q, and Ax, Av, At. The results of a few of these runs are pre- 
sented on the following pages. We normalize the velocity domain by choosing Q — 1 
so that velocities are only allowed in the interval [—1, 1]. In each of the runs, it is 
important to properly balance the choices of T and L. If L is taken too small or T 
too large, the computations of the field at each gridpoint will be invalid after some 
small time To < T. Thus, we must ensure that the domain of validity contains a 
region [0, T] x [—/,/] for some I < L (see Figure [T]). For each of the following runs, 
the mesh sizes are taken to be At = Ax = Av = 0.01 and L = 50 so that the 
computations are valid over a reasonably long time interval. More specifically, since 
maximum velocities are on average near Q — 1, we expect the time of validity to 



displays the computed electric field and 
in the next section will display the same 



be approximately T = L/Q w 50. Figure 2 
conserved energy for T = 30, while Figure 3 
information for T = 50. 

From Figure |2j we can see that the computed electric field is oscillatory, but with 
decreasing amplitude and mean for < t < 30. In addition, after considering values 
of t near which the peaks of these oscillations occur, we may conclude this amplitude 
decreases around a rate of t~ l for t e [0, 30]. This result is demonstrated in Table [2j 
Specifically, the last row of the table leads us to believe that for t G [0, 30] , 

sup|£(t,^)| « 1.1 x 10~ 4 r 1 . 

To the right of the computed electric field, Figure [2] shows the net energy in 
the computational domain for times < t < 30. The deviation of the computed 
energy from its initially computed value is typically seen as a good measure of the 
accuracy of the method over time. In this case, the values of the energy range between 
1.747 x 10 -2 , initially, and 1.665 x 10 -2 , near time T = 30. The relative change is 
calculated as 

1.755 x Hr 2 - 1.657 x 10~ 2 

%change = « 4.92%. 

1.657 x 10 1 

Thus, we expect our computations to be within 4% and 5% of their actual value. 
Breakdown Outside the Domain of Validity 

While the particle method does yield accurate and efficient approximations of 0, 
limitations do exist within the formulation and implementation. Of course, a limita- 
tion of any approximation of a particle distribution without compact support is the 
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5 10 15 20 25 30 35 40 45 50 



Figure 3: The computed electric field (left, lined), asymptotic behavior (left, dashed), and net energy 
(right) for < t < 50 

truncation of the spatial and velocity domains. One can only simulate domains of 
finite length, and the previously described particle method will determine the true 
effects observed, but only within the domain of validity. However, the size of the 
domain of validity is another limitation of the method. In Figure [3j we can see that 
after time t = 30, the behavior of the calculated field changes slightly. The decay to 
zero is impeded and the mean of the oscillation begins to increase. In addition, the 
figure demonstrates that the energy begins to decrease further away from its initial 
value after time 30, as well. Within the simulation, particle velocities have increased, 
causing the expected domain of validity to shrink. Thus, the computations are only 
valid up to a time which is strictly less than the stopping time, after which other par- 
ticles beginning outside of the truncated spatial domain will influence field behavior 
within. This displays the limitation of the particle method, graphically represented 
in the previous section by Figure [TJ Computations will inevitably fail to be accurate 
on any portion of the spatial domain after some time. The only remedy within the 
context of a particle method is to enlarge the initial spatial domain, resulting in more 
expensive computations. 
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